classdef cASL < cMRI
    % CASL The class cASL is an object of ASL data
    %   Its properties include:
    %       dASL    .....   ASL data (substraction)
    %       dMean   .....   The averaged ASL
    %       dM0     .....   The first repetition (used for M0)
    %       vLabel  .....   The start repetition of labeled MRI
    %       vControl.....   The start repetition of control MRI
    %       TI2     .....   The capturing time
    %       Tau     .....   The bolus length
    %   NOTE: There are other properties reminded to be added in future,
    %   such as labeling efficiency, Tissue spinlattice relaxation
    %   constant, and so on.
    %
    %   Its methods include:
    %       cASL(varargin)  ... Constructor
    %       setM0(obj)      ... Compute M0 then save it with cMRI format
    %       setASL(obj)     ... Compute ASL then save it with cMRI format
    %       getASL(obj,varargin)... Get data or dim of dASL
    %       getMean(obj,varargin).. Get data or dim of dMean.
    %   NOTE: You might be interested in adding your own methods...
    %    
    % COPYRIGHT RESERVED 2012 @ INRIA/IRISA VISAGE, LEI YU
    % See also: cMRI, cDIVASL
    %
    
    properties
        % Data that saved
        dASL = [];  % Asl Data
        dMean = []; % Averaged Asl data
        dM0 = [];
        dGray = []; % Template of gray matter
        dWhite = []; % Template of white matter
        dShell = []; % Template of brain shell
        
        dCBF = cMRI; % The CBF
        
        % Some useful values
        vLabel = 3; % Start repetition of Label image
        vControl = 2; % Start repetition of Control image
        vM0 = 1; % The repetition of M0
        
        nPairs = 0;
        
        % Model parameters, varying with different machines
        pAlpha = 0.95;
        pLambda = 0.9; % ml.g(-1)
        pFactor = 1000*3; % this value depends on the units
        pTI2 = 0;    % Should be set manually according the sampling time
        pT1b = 1.5;
        pTau = 0.7; % in seconds
    end
    
    methods
        function obj = cASL(varargin)
            % Construction function
            % the input should be a cell or array of cell
            
            if nargin ~=0 % Provide default settings
                if length(varargin) == 1
                    varargin = varargin{:};
                    inputs = varargin;
                    obj(length(inputs),1) = cASL;
                    for i = 1:length(inputs)
                        obj(i) = obj(i).initialMRI(inputs{i});
                        if nargin >= 2
                            obj(i).setParameters(varargin{2:end});
                        end
                        obj(i) = obj(i).setM0();
                        obj(i) = obj(i).setASL();
                    end
                elseif length(varargin) >= 2 %% if we have m0 and asl
                    f_list_m0 = varargin{1};
                    f_list_asl = varargin{2};
                    f_list_mri= varargin{3};
                    obj(length(f_list_m0),1) = cASL;
                    for i = 1:length(f_list_m0)
                        obj(i) = obj(i).initialMRI(f_list_mri{i});
                        obj(i) = obj(i).setM0(f_list_m0{i});
                        obj(i) = obj(i).setASL(f_list_asl{i});
                    end
                    
                end
            else
            end
            
        end
    end
    
    methods (Access = private)
        % This function is to set the M0 image function (from data or file)
        function obj = setM0(obj,varargin)
            if isempty(varargin)
                obj.dM0 = cMRI(obj.dMRI(:,:,:,obj.vM0));
            else
                obj.dM0 = cMRI(varargin{1});
            end
        end
        % This function is to set the ASL image function (from data or file)
        function obj = setASL(obj,varargin)
            if isempty(varargin)
                % Compute ASL data (substraction) and set the ASL infos.
                idLabel = obj.vLabel:2:obj.rImg;
                idControl = obj.vControl:2:obj.rImg;
                if ~isequal(length(idLabel),length(idControl))
                    error('The volumns of label and control images are not identical.');
                end
                obj.nPairs = length(idLabel);
                obj.dASL = cMRI(obj.dMRI(:,:,:,idLabel) - obj.dMRI(:,:,:,idControl));
%                 obj.dMean = cMRI(sum(obj.dASL.getImg('MRI'),4)/obj.nPairs);
            else
                obj.dASL = cMRI(varargin{1});
                obj.nPairs = size(obj.dASL.getImg('MRI'),4);
            end
            obj.dMean = cMRI(sum(obj.dASL.getImg('MRI'),4)/obj.nPairs);
        end
        
    end
    methods
        % This function is to get asl image with single repetition
        function varargout = getASL(obj, varargin)
            varargout{1} = obj.dASL.getImg(varargin);
            
        end
        
        % This function is to get mean value of the 
        function varargout = getMean(obj, varargin)
            switch lower(varargin{1})
                case 'def'
                    varargout{1} = mean(obj.dASL.getImg('mri',varargin{2},varargin{3},varargin{4},varargin{5}),4);
                case 'mri'
                varargout{1} = obj.dMean.getImg(varargin);
            end
        end
        
        function varargout = getM0(obj, varargin)
            varargout{1} = obj.dM0.getImg(varargin);
        end
        
        % This function is to set the Template
        function obj = setTemplate(obj, varargin)
            varargin = varargin{:};
            if ~isempty(varargin)
                nVar = length(varargin);
                i = 1;
                while i<=nVar-1
                    switch lower(varargin{i})
                        case 'gray'
                            i = i + 1;
                            if ischar(varargin{i})
                                obj.dGray = cMRI(varargin{i});i = i + 1;
                            else
                                obj.dGray = varargin{i};i = i+ 1;
                            end
                        case 'white'
                            i = i + 1;
                            if ischar(varargin{i})
                                obj.dWhite = cMRI(varargin{i});i = i + 1;
                            else
                                obj.dWhite = varargin{i};i = i+ 1;
                            end
                        case 'shell'
                            i = i + 1;
                            if ischar(varargin{i})
                                obj.dShell = cMRI(varargin{i});i = i + 1;
                            else
                                obj.dShell = varargin{i};i = i+ 1;
                            end
                        otherwise
                            error('Failure');
                    end
                    
                end
            else
                error( ' Template setting failure');
            end
        end
        
        % This function is to get the masked ASL or Mean or M0
        function varargout = getMaskedImg(obj,dataType,maskType, varargin)
            switch lower(maskType)
                case 'gray'
                    maskOut = obj.dGray;
                case 'white'
                    maskOut = obj.dWhite;
                case 'shell'
                    maskOut = obj.dShell;
                otherwise
                    error('Incorrect mask type');
            end
            
            switch lower(dataType)
                case 'asl'
                    dataOut = obj.dASL;
                case 'mean'
                    dataOut = obj.dMean;
                case 'm0'
                    dataOut = obj.dM0;
                otherwise
                    error('Incorrect data type');
            end
            
            
            % If there is slice or repetions specified
            if ~isempty(varargin)
%                 varargin = varargin{:};
                slices = varargin{1};
                rep = varargin{2};
                data2Mask = dataOut.getImg('mri',:,:,slices,rep);
                dataMask = maskOut.getImg('mri',:,:,slices,rep)>0.5;
            else
                data2Mask = dataOut.getImg('mri');
                dataMask = maskOut.getImg('mri')>0.5;
            end
            
%             dataMask = maskOut.getImg('mri') > 0.5;
            
            
            [w,h,slice,rep] = size(data2Mask);
            dataMasked = zeros(w,h,slice,rep);
            for i = 1:rep
                dataMasked(:,:,:,i) = data2Mask(:,:,:,i).*dataMask;
            end
            
%             dataMasked = obj.maskImg(data2Mask,maskOut);
            
            varargout{1} = dataMasked;
            varargout{2} = dataMask;
            varargout{3} = data2Mask;
                
        end
        
        % This function is to calculate the CBF using single TI with
        % buxton's model
        function obj = calCBFImg(obj,varargin)
            
            i = 1;commands = [];
            while i<length(varargin)
                if ischar(varargin{i})
                    commands = [commands,varargin{i}];
                end
                i = i+ 1;
            end
            
            if isempty(commands)
                commands = {'buxton','mask','gray'};
            end
            
           
            
            
            isMask = cell2mat(regexpi(commands,{'mask'}));
            if isMask
                maskType = (regexpi([commands{:}],{'gray|shell|white'},'match'));
                M0 = obj.getMaskedImg('m0',maskType{1}{1});
                Mean = obj.getMaskedImg('mean',maskType{1}{1});
            else
                M0 = obj.getM0('mri');
                Mean = obj.getMean('mri');
            end
              
            isBuxton = cell2mat(regexpi(commands,{'buxton'}));
            if isBuxton
                slices = size(Mean,3);
                CBF = zeros(size(Mean));
                for i = 1:slices
                    CBF(:,:,i) = obj.calBuxtonCBF(Mean(:,:,i),M0(:,:,i),i);
                end
                obj.dCBF = obj.dCBF.initialMRI(CBF);
            end
            
        end
        
        % This function is to calculate the CBF for one slice of image
        % So the input data dataIn should be 64x64
        function CBF = calBuxtonCBF(obj,dataMean,dataM0,pSlice)
%             CBF = zeros(size(dataMean));
            factors = obj.pFactor*...
                      obj.pLambda/...
                      obj.pAlpha/...
                      obj.pTau/...
                      exp(-(obj.pTI2+(pSlice-1)*0.045)/obj.pT1b);
                  
            CBF = factors*dataMean./(dataM0+eps); % to prevent NAN values
        end
        
    end

end

